function [ omegaENNRate ] = transraterate( vN, vNRate, L, h, RE, eSq )
%TRANSRATERATE 求解位移角速度的变化率
%  
% Tips:
% 1. 导航系N使用ENU地理坐标系
%
% Input Arguments:
% # vN: 3*1向量，N系下的地速，单位m/s
% # vNRate:  3*1向量，N系下的加速度，单位m/s
% # L: 标量，纬度，单位rad
% # h: 标量，高度，单位m
% # RE: 标量，地球长半轴，单位m
% # eSq: 标量，地球偏心率的平方，无单位
%
% Output Arguments:
% # omegaENNRate：3*1向量，位移角速度的变化率，单位rad/s/s
% 
% References:
% # 算法文档-轨迹发生器输出——辅助输出

[RM, RN] = earthradiiofcurv(L, RE, eSq);
[LRate, lambdaRate] = latlonrate(vN(1, 1), vN(2, 1), L, h, RE, eSq);

% 子午圈、卯酉圈曲率半径随纬度的变换率 
k = 1 - eSq * (sin(L))^2;
RMRate = 3 * RE * eSq * (1-eSq) * sin(L) * cos(L) * LRate / k^(5/2); % Eq. 3.1.1
RNRate = RE * eSq * sin(L) * cos(L) * LRate / k^(3/2);

% 纬度、经度变化率的微分
LRateRate = vNRate(2, 1) / (RM + h) - vN(2, 1)*(RMRate + vN(3, 1)) / (RM + h)^2; % Eq. 21.3.19
lambdaRateRate = vNRate(1, 1) / (RN + h) / cos(L) - vN(1, 1) * ((RNRate+vN(3, 1))*cos(L) - (RN + h)*sin(L) * LRate) / ((RN + h) * cos(L))^2;

omegaENNRate = [-LRateRate, lambdaRateRate*cos(L)-lambdaRate*sin(L)*LRate,  lambdaRateRate*sin(L)+lambdaRate*cos(L)*LRate]';
end

